# Set the working direction to where your files are located.
setwd('C:\\Users\\mawal\\OneDrive - Binghamton University\\Desktop\\Desktop_Folders\\Upwork\\Kreps_2')
df = read.csv('study_2_.csv')[-1] # File manipulated in python for study 2

library(dplyr) # Needed for data manipulation

# Changing the names of the treatments to make them more intuitive
df = df %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_12', 'hum_climate', FL_4_DO)) %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_13', 'hum_economy', treatments)) %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_14', 'hum_business', treatments)) %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_15', 'gpt_climate', treatments)) %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_16', 'gpt_economy', treatments)) %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_17', 'gpt_business', treatments))

# Turning the treaments variable into a factor variable
df$treatments = as.factor(df$treatments)

# Turning Integer into numeric
df$Q384_4 = as.numeric(df$Q384_4)

###########################################################
# Plots with Q384_4 as the DV and treaments as the IV
###########################################################
# Running the regression with the Slider as the DV
regr = lm(Q384_4 ~ treatments, data = df)
summary(regr)

library(MASS) # to use the mvrnorm function

# Creating a multivariate normal distribution with mu and sigma set from the regression
coefs = mvrnorm(n = 10000, mu = coefficients(regr), Sigma = vcov(regr))

# Doing some diagnositcs by comparing the newly created data with the regression output
coefficients(regr) #coefficients
colMeans(coefs) #Means of the columns to match coefficients

sqrt(diag(vcov(regr))) #errors
apply(coefs, 2, sd) # errors to match errors from the regression

#Turning coefs into data frame
coefs = as.data.frame(coefs)

test1 = list(1) # Temporary place holder
test1 = as.data.frame(test1)
names(test1)[names(test1) == "X1"] <- "lower" # 5th percentile placeholder
test1$median = 0 # 50th percentile placeholder
test1$upper = 0 # 95th percentile placeholder

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 samples
  coefs$xb[row] = coefs$`(Intercept)`[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test1 = rbind(test1, vec1) # Binding the new row to the main test data frame


test1 = test1[-1,] # Remove the first row, just placeholder data
test1$type = 'GPT3' # Determining the type of the base category
test1$topic = 'Local Business' # Dertermining the topic of the base category


test2 = list(1)
test2 = as.data.frame(test2)
names(test2)[names(test2) == "X1"] <- "lower"
test2$median = 0
test2$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test2 = rbind(test2, vec1) # Binding the new row to the main test data frame

test2 = test2[-1,]
test2$type = 'GPT3'
test2$topic = 'Climate'

test3 = list(1)
test3 = as.data.frame(test3)
names(test3)[names(test3) == "X1"] <- "lower"
test3$median = 0
test3$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test3 = rbind(test3, vec1) # Binding the new row to the main test data frame

test3 = test3[-1,]
test3$type = 'GPT3'
test3$topic = 'Economy'


test4 = list(1)
test4 = as.data.frame(test4)
names(test4)[names(test4) == "X1"] <- "lower"
test4$median = 0
test4$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_business[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test4 = rbind(test4, vec1) # Binding the new row to the main test data frame

test4 = test4[-1,]
test4$type = 'Human'
test4$topic = 'Local Business'


test5 = list(1)
test5 = as.data.frame(test5)
names(test5)[names(test5) == "X1"] <- "lower"
test5$median = 0
test5$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test5 = rbind(test5, vec1) # Binding the new row to the main test data frame

test5 = test5[-1,]
test5$type = 'Human'
test5$topic = 'Climate'


test6 = list(1)
test6 = as.data.frame(test6)
names(test6)[names(test6) == "X1"] <- "lower"
test6$median = 0
test6$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test6 = rbind(test6, vec1) # Binding the new row to the main test data frame

test6 = test6[-1,]
test6$type = 'Human'
test6$topic = 'Economy'

# Combining all of the simulation into one dataframe
final = rbind(test1, test2, test3, test4, test5, test6)

library(ggplot2) # Necessary to create the plot
# Creating the plot 
p<- ggplot(final, aes(x=topic, y=median, col = type )) + 
  geom_point(position=position_dodge(0.25))+
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2,
                position=position_dodge(0.25))
slider_plot = p+labs(title="Spending Money on Environment/Economy", x="Treatment Group", y = "Amount out of $100")+
  theme_classic() +
  scale_color_manual(values=c('black','gray')) +
  labs(col = "GPT3 vs Human")
slider_plot

# Saving the plot to your working directory
ggsave('slider_plot.pdf', slider_plot)

##########################################################
###########################################################
# Plots with political behavior count as the DV and treaments as the IV
###########################################################
rm(list=ls()) # Wiping out the global environment
setwd('C:\\Users\\mawal\\OneDrive - Binghamton University\\Desktop\\Desktop_Folders\\Upwork\\Kreps_2')
df = read.csv('data_2.csv')[-1] # Contains the political behavior variables

# Changing the names of the treatments to make them more intuitive
library(dplyr)
df = df %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_12', 'hum_climate', FL_4_DO)) %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_13', 'hum_economy', treatments)) %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_14', 'hum_business', treatments)) %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_15', 'gpt_climate', treatments)) %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_16', 'gpt_economy', treatments)) %>%
  mutate(treatments = ifelse(FL_4_DO == 'FL_17', 'gpt_business', treatments))

# Turning political behavior into a numeric variable
df$pol_behaviour = as.numeric(df$pol_behaviour)

# Running regressions with behavior count as the DV and treatments as the IV
regr = lm(behavior_count ~ treatments, data = df)
summary(regr)

library(MASS)
#Using the variance-covariance matrix and coefficients to create 10000 multivariate normal distributions.
coefs = mvrnorm(n = 10000, mu = coefficients(regr), Sigma = vcov(regr))

# Doing some diagnositcs by comparing the newly created data with the regression output
coefficients(regr) #coefficients
colMeans(coefs)

sqrt(diag(vcov(regr))) #errors
apply(coefs, 2, sd)

#Turning coefs into data frame
coefs = as.data.frame(coefs)

test1 = list(1)
test1 = as.data.frame(test1)
names(test1)[names(test1) == "X1"] <- "lower"
test1$median = 0
test1$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test1 = rbind(test1, vec1) # Binding the new row to the main test data frame


test1 = test1[-1,]
test1$type = 'Human'
test1$topic = 'Local Business'


test2 = list(1)
test2 = as.data.frame(test2)
names(test2)[names(test2) == "X1"] <- "lower"
test2$median = 0
test2$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test2 = rbind(test2, vec1) # Binding the new row to the main test data frame

test2 = test2[-1,]
test2$type = 'GPT3'
test2$topic = 'Climate'


test3 = list(1)
test3 = as.data.frame(test3)
names(test3)[names(test3) == "X1"] <- "lower"
test3$median = 0
test3$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test3 = rbind(test3, vec1) # Binding the new row to the main test data frame


test3 = test3[-1,]
test3$type = 'GPT3'
test3$topic = 'Economy'


test4 = list(1)
test4 = as.data.frame(test4)
names(test4)[names(test4) == "X1"] <- "lower"
test4$median = 0
test4$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_business[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test4 = rbind(test4, vec1) # Binding the new row to the main test data frame

test4 = test4[-1,]
test4$type = 'GPT3'
test4$topic = 'Local Business'


test5 = list(1)
test5 = as.data.frame(test5)
names(test5)[names(test5) == "X1"] <- "lower"
test5$median = 0
test5$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test5 = rbind(test5, vec1) # Binding the new row to the main test data frame


test5 = test5[-1,]
test5$type = 'Human'
test5$topic = 'Climate'


test6 = list(1)
test6 = as.data.frame(test6)
names(test6)[names(test6) == "X1"] <- "lower"
test6$median = 0
test6$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test6 = rbind(test6, vec1) # Binding the new row to the main test data frame

test6 = test6[-1,]
test6$type = 'Human'
test6$topic = 'Economy'

# Binding the dataframes together
final = rbind(test1, test2, test3, test4, test5, test6)

# Plotting the results
library(ggplot2)
p<- ggplot(final, aes(x=topic, y=median, col = type )) + 
  geom_point(position=position_dodge(0.25))+
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2,
                position=position_dodge(0.25))
behavior_regr_plot = p+labs(title="Involvement in Politics", x="Treatment Group", y = "Number of Activities Scale")+
  theme_classic() +
  scale_color_manual(values=c('black','gray')) +
  labs(col = "GPT3 vs Human")
behavior_regr_plot

ggsave('behavior_regr_plot.pdf', behavior_regr_plot)


#################################################
###########################################################
# Plots with Attend the DV and treaments as the IV
###########################################################
# Running the regression
regr = glm(attend ~ treatments, data = df, family = 'binomial')
summary(regr)

library(MASS)
#Using the variance-covariance matrix and coefficients to create 10000 multivariate normal distributions.
coefs = mvrnorm(n = 10000, mu = coefficients(regr), Sigma = vcov(regr))

# Doing some diagnositcs by comparing the newly created data with the regression output
coefficients(regr) #coefficients
colMeans(coefs)

sqrt(diag(vcov(regr))) #errors
apply(coefs, 2, sd)

#Turning coefs into data frame
coefs = as.data.frame(coefs)

test1 = list(1)
test1 = as.data.frame(test1)
names(test1)[names(test1) == "X1"] <- "lower"
test1$median = 0
test1$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test1 = rbind(test1, vec1) # Binding the new row to the main test data frame


test1 = test1[-1,]
test1$type = 'Human'
test1$topic = 'Local Business'

test1$p = 1/(1+exp(-(test1$median))) # Using the logit equation to turn into probabilities.
test1$pub = 1/(1+ exp(-(test1$upper))) # Using the logit equation to turn into probabilities.
test1$plb = 1/(1+ exp(-(test1$lower))) # Using the logit equation to turn into probabilities.


test2 = list(1)
test2 = as.data.frame(test2)
names(test2)[names(test2) == "X1"] <- "lower"
test2$median = 0
test2$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test2 = rbind(test2, vec1) # Binding the new row to the main test data frame

test2 = test2[-1,]

test2$type = 'GPT3'
test2$topic = 'Climate'

test2$p = 1/(1+exp(-(test2$median)))
test2$pub = 1/(1+ exp(-(test2$upper)))
test2$plb = 1/(1+ exp(-(test2$lower)))


test3 = list(1)
test3 = as.data.frame(test3)
names(test3)[names(test3) == "X1"] <- "lower"
test3$median = 0
test3$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test3 = rbind(test3, vec1) # Binding the new row to the main test data frame


test3 = test3[-1,]
test3$type = 'GPT3'
test3$topic = 'Economy'

test3$p = 1/(1+exp(-(test3$median)))
test3$pub = 1/(1+ exp(-(test3$upper)))
test3$plb = 1/(1+ exp(-(test3$lower)))


test4 = list(1)
test4 = as.data.frame(test4)
names(test4)[names(test4) == "X1"] <- "lower"
test4$median = 0
test4$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_business[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test4 = rbind(test4, vec1) # Binding the new row to the main test data frame

test4 = test4[-1,]
test4$type = 'GPT3'
test4$topic = 'Local Business'

test4$p = 1/(1+exp(-(test4$median)))
test4$pub = 1/(1+ exp(-(test4$upper)))
test4$plb = 1/(1+ exp(-(test4$lower)))


test5 = list(1)
test5 = as.data.frame(test5)
names(test5)[names(test5) == "X1"] <- "lower"
test5$median = 0
test5$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test5 = rbind(test5, vec1) # Binding the new row to the main test data frame


test5 = test5[-1,]
test5$type = 'Human'
test5$topic = 'Climate'

test5$p = 1/(1+exp(-(test5$median)))
test5$pub = 1/(1+ exp(-(test5$upper)))
test5$plb = 1/(1+ exp(-(test5$lower)))


test6 = list(1)
test6 = as.data.frame(test6)
names(test6)[names(test6) == "X1"] <- "lower"
test6$median = 0
test6$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test6 = rbind(test6, vec1) # Binding the new row to the main test data frame

test6 = test6[-1,]
test6$type = 'Human'
test6$topic = 'Economy'

test6$p = 1/(1+exp(-(test6$median)))
test6$pub = 1/(1+ exp(-(test6$upper)))
test6$plb = 1/(1+ exp(-(test6$lower)))


final = rbind(test1, test2, test3, test4, test5, test6)

library(ggplot2)
p<- ggplot(final, aes(x=topic, y=p, col = type )) + 
  geom_point(position=position_dodge(0.25))+
  geom_errorbar(aes(ymin=plb, ymax=pub), width=.2,
                position=position_dodge(0.25))
attend_plot = p+labs(title="Local Government Attendance", x="Treatment Group", y = "Probability of Attending")+
  theme_classic() +
  scale_color_manual(values=c('black','gray')) +
  labs(col = "GPT3 vs Human")
attend_plot

ggsave('attend_plot.pdf', attend_plot)
###################################################################
###########################################################
# Plots with Participate as the DV and treaments as the IV
###########################################################
# Regressions
regr = glm(parti ~ treatments, data = df, family = 'binomial')
summary(regr)

library(MASS)
#Using the variance-covariance matrix and coefficients to create 10000 multivariate normal distributions.
coefs = mvrnorm(n = 10000, mu = coefficients(regr), Sigma = vcov(regr))

# Doing some diagnositcs by comparing the newly created data with the regression output
coefficients(regr) #coefficients
colMeans(coefs)

sqrt(diag(vcov(regr))) #errors
apply(coefs, 2, sd)

#Turning coefs into data frame
coefs = as.data.frame(coefs)

test1 = list(1)
test1 = as.data.frame(test1)
names(test1)[names(test1) == "X1"] <- "lower"
test1$median = 0
test1$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test1 = rbind(test1, vec1) # Binding the new row to the main test data frame


test1 = test1[-1,]
test1$type = 'Human'
test1$topic = 'Local Business'

test1$p = 1/(1+exp(-(test1$median)))
test1$pub = 1/(1+ exp(-(test1$upper)))
test1$plb = 1/(1+ exp(-(test1$lower)))


test2 = list(1)
test2 = as.data.frame(test2)
names(test2)[names(test2) == "X1"] <- "lower"
test2$median = 0
test2$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test2 = rbind(test2, vec1) # Binding the new row to the main test data frame

test2 = test2[-1,]
test2$type = 'GPT3'
test2$topic = 'Climate'

test2$p = 1/(1+exp(-(test2$median)))
test2$pub = 1/(1+ exp(-(test2$upper)))
test2$plb = 1/(1+ exp(-(test2$lower)))


test3 = list(1)
test3 = as.data.frame(test3)
names(test3)[names(test3) == "X1"] <- "lower"
test3$median = 0
test3$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test3 = rbind(test3, vec1) # Binding the new row to the main test data frame


test3 = test3[-1,]
test3$type = 'GPT3'
test3$topic = 'Economy'

test3$p = 1/(1+exp(-(test3$median)))
test3$pub = 1/(1+ exp(-(test3$upper)))
test3$plb = 1/(1+ exp(-(test3$lower)))


test4 = list(1)
test4 = as.data.frame(test4)
names(test4)[names(test4) == "X1"] <- "lower"
test4$median = 0
test4$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_business[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test4 = rbind(test4, vec1) # Binding the new row to the main test data frame

test4 = test4[-1,]
test4$type = 'GPT3'
test4$topic = 'Local Business'

test4$p = 1/(1+exp(-(test4$median)))
test4$pub = 1/(1+ exp(-(test4$upper)))
test4$plb = 1/(1+ exp(-(test4$lower)))



test5 = list(1)
test5 = as.data.frame(test5)
names(test5)[names(test5) == "X1"] <- "lower"
test5$median = 0
test5$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test5 = rbind(test5, vec1) # Binding the new row to the main test data frame


test5 = test5[-1,]
test5$type = 'Human'
test5$topic = 'Climate'

test5$p = 1/(1+exp(-(test5$median)))
test5$pub = 1/(1+ exp(-(test5$upper)))
test5$plb = 1/(1+ exp(-(test5$lower)))


test6 = list(1)
test6 = as.data.frame(test6)
names(test6)[names(test6) == "X1"] <- "lower"
test6$median = 0
test6$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test6 = rbind(test6, vec1) # Binding the new row to the main test data frame

test6 = test6[-1,]
test6$type = 'Human'
test6$topic = 'Economy'

test6$p = 1/(1+exp(-(test6$median)))
test6$pub = 1/(1+ exp(-(test6$upper)))
test6$plb = 1/(1+ exp(-(test6$lower)))


final = rbind(test1, test2, test3, test4, test5, test6)

library(ggplot2)
p<- ggplot(final, aes(x=topic, y=p, col = type )) + 
  geom_point(position=position_dodge(0.25))+
  geom_errorbar(aes(ymin=plb, ymax=pub), width=.2,
                position=position_dodge(0.25))
participation_plot = p+labs(title="Participation in Political Protest", x="Treatment Group", y = "Probability of Participating")+
  theme_classic() +
  scale_color_manual(values=c('black','gray')) +
  labs(col = "GPT3 vs Human")
participation_plot

ggsave('participation_plot.pdf', participation_plot)
###################################################################
###########################################################
# Plots with Donate the DV and treaments as the IV
###########################################################
#Regressions
regr = glm(donate ~ treatments, data = df, family = 'binomial')
summary(regr)

library(MASS)
#Using the variance-covariance matrix and coefficients to create 10000 multivariate normal distributions.
coefs = mvrnorm(n = 10000, mu = coefficients(regr), Sigma = vcov(regr))

# Doing some diagnositcs by comparing the newly created data with the regression output
coefficients(regr) #coefficients
colMeans(coefs)

sqrt(diag(vcov(regr))) #errors
apply(coefs, 2, sd)

#Turning coefs into data frame
coefs = as.data.frame(coefs)

test1 = list(1)
test1 = as.data.frame(test1)
names(test1)[names(test1) == "X1"] <- "lower"
test1$median = 0
test1$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test1 = rbind(test1, vec1) # Binding the new row to the main test data frame


test1 = test1[-1,]
test1$type = 'GPT3'
test1$topic = 'Economy'

test1$p = 1/(1+exp(-(test1$median)))
test1$pub = 1/(1+ exp(-(test1$upper)))
test1$plb = 1/(1+ exp(-(test1$lower)))


test2 = list(1)
test2 = as.data.frame(test2)
names(test2)[names(test2) == "X1"] <- "lower"
test2$median = 0
test2$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test2 = rbind(test2, vec1) # Binding the new row to the main test data frame

test2 = test2[-1,]
test2$type = 'GPT3'
test2$topic = 'Climate'

test2$p = 1/(1+exp(-(test2$median)))
test2$pub = 1/(1+ exp(-(test2$upper)))
test2$plb = 1/(1+ exp(-(test2$lower)))



test3 = list(1)
test3 = as.data.frame(test3)
names(test3)[names(test3) == "X1"] <- "lower"
test3$median = 0
test3$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_business[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test3 = rbind(test3, vec1) # Binding the new row to the main test data frame


test3 = test3[-1,]
test3$type = 'Human'
test3$topic = 'Local Business'

test3$p = 1/(1+exp(-(test3$median)))
test3$pub = 1/(1+ exp(-(test3$upper)))
test3$plb = 1/(1+ exp(-(test3$lower)))


test4 = list(1)
test4 = as.data.frame(test4)
names(test4)[names(test4) == "X1"] <- "lower"
test4$median = 0
test4$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_business[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test4 = rbind(test4, vec1) # Binding the new row to the main test data frame

test4 = test4[-1,]
test4$type = 'GPT3'
test4$topic = 'Local Business'

test4$p = 1/(1+exp(-(test4$median)))
test4$pub = 1/(1+ exp(-(test4$upper)))
test4$plb = 1/(1+ exp(-(test4$lower)))


test5 = list(1)
test5 = as.data.frame(test5)
names(test5)[names(test5) == "X1"] <- "lower"
test5$median = 0
test5$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test5 = rbind(test5, vec1) # Binding the new row to the main test data frame


test5 = test5[-1,]
test5$type = 'Human'
test5$topic = 'Climate'

test5$p = 1/(1+exp(-(test5$median)))
test5$pub = 1/(1+ exp(-(test5$upper)))
test5$plb = 1/(1+ exp(-(test5$lower)))


test6 = list(1)
test6 = as.data.frame(test6)
names(test6)[names(test6) == "X1"] <- "lower"
test6$median = 0
test6$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test6 = rbind(test6, vec1) # Binding the new row to the main test data frame

test6 = test6[-1,]
test6$type = 'Human'
test6$topic = 'Economy'

test6$p = 1/(1+exp(-(test6$median)))
test6$pub = 1/(1+ exp(-(test6$upper)))
test6$plb = 1/(1+ exp(-(test6$lower)))


final = rbind(test1, test2, test3, test4, test5, test6)

library(ggplot2)
p<- ggplot(final, aes(x=topic, y=p, col = type )) + 
  geom_point(position=position_dodge(0.25))+
  geom_errorbar(aes(ymin=plb, ymax=pub), width=.2,
                position=position_dodge(0.25))
donate_plot = p+labs(title="Donate to Political Campaign", x="Treatment Group", y = "Probability of Donating")+
  theme_classic() +
  scale_color_manual(values=c('black','gray')) +
  labs(col = "GPT3 vs Human")
donate_plot

ggsave('donate_plot.pdf', donate_plot)
###################################################################
###########################################################
# Plots with Write as the DV and treaments as the IV
###########################################################

regr = glm(write ~ treatments, data = df, family = 'binomial')
summary(regr)

library(MASS)
#Using the variance-covariance matrix and coefficients to create 10000 multivariate normal distributions.
coefs = mvrnorm(n = 10000, mu = coefficients(regr), Sigma = vcov(regr))

# Doing some diagnositcs by comparing the newly created data with the regression output
coefficients(regr) #coefficients
colMeans(coefs)

sqrt(diag(vcov(regr))) #errors
apply(coefs, 2, sd)

#Turning coefs into data frame
coefs = as.data.frame(coefs)

test1 = list(1)
test1 = as.data.frame(test1)
names(test1)[names(test1) == "X1"] <- "lower"
test1$median = 0
test1$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test1 = rbind(test1, vec1) # Binding the new row to the main test data frame


test1 = test1[-1,]
test1$type = 'GPT3'
test1$topic = 'Climate'

test1$p = 1/(1+exp(-(test1$median)))
test1$pub = 1/(1+ exp(-(test1$upper)))
test1$plb = 1/(1+ exp(-(test1$lower)))


test2 = list(1)
test2 = as.data.frame(test2)
names(test2)[names(test2) == "X1"] <- "lower"
test2$median = 0
test2$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test2 = rbind(test2, vec1) # Binding the new row to the main test data frame

test2 = test2[-1,]

test2$type = 'GPT3'
test2$topic = 'Economy'

test2$p = 1/(1+exp(-(test2$median)))
test2$pub = 1/(1+ exp(-(test2$upper)))
test2$plb = 1/(1+ exp(-(test2$lower)))


test3 = list(1)
test3 = as.data.frame(test3)
names(test3)[names(test3) == "X1"] <- "lower"
test3$median = 0
test3$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_business[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test3 = rbind(test3, vec1) # Binding the new row to the main test data frame


test3 = test3[-1,]
test3$type = 'Human'
test3$topic = 'Local Business'

test3$p = 1/(1+exp(-(test3$median)))
test3$pub = 1/(1+ exp(-(test3$upper)))
test3$plb = 1/(1+ exp(-(test3$lower)))


test4 = list(1)
test4 = as.data.frame(test4)
names(test4)[names(test4) == "X1"] <- "lower"
test4$median = 0
test4$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentsgpt_business[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test4 = rbind(test4, vec1) # Binding the new row to the main test data frame

test4 = test4[-1,]
test4$type = 'GPT3'
test4$topic = 'Local Business'

test4$p = 1/(1+exp(-(test4$median)))
test4$pub = 1/(1+ exp(-(test4$upper)))
test4$plb = 1/(1+ exp(-(test4$lower)))


test5 = list(1)
test5 = as.data.frame(test5)
names(test5)[names(test5) == "X1"] <- "lower"
test5$median = 0
test5$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_climate[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test5 = rbind(test5, vec1) # Binding the new row to the main test data frame


test5 = test5[-1,]
test5$type = 'Human'
test5$topic = 'Climate'

test5$p = 1/(1+exp(-(test5$median)))
test5$pub = 1/(1+ exp(-(test5$upper)))
test5$plb = 1/(1+ exp(-(test5$lower)))


test6 = list(1)
test6 = as.data.frame(test6)
names(test6)[names(test6) == "X1"] <- "lower"
test6$median = 0
test6$upper = 0

for(row in 1:10000) { # We want to take the median and qi's from 10,0000 sammples
  coefs$xb[row] = coefs$`(Intercept)`[row] + coefs$treatmentshum_economy[row]
}
qi = coefs %>% #T  Taking quantities of interest from the xb variable
  summarize(quants = quantile(xb, probs = c(0.05, .5, 0.95)))
vec1 = c(qi[1,1], qi[2,1], qi[3,1]) #Turning these quantities into a list
test6 = rbind(test6, vec1) # Binding the new row to the main test data frame

test6 = test6[-1,]
test6$type = 'Human'
test6$topic = 'Economy'

test6$p = 1/(1+exp(-(test6$median)))
test6$pub = 1/(1+ exp(-(test6$upper)))
test6$plb = 1/(1+ exp(-(test6$lower)))


final = rbind(test1, test2, test3, test4, test5, test6)

library(ggplot2)
p<- ggplot(final, aes(x=topic, y=p, col = type )) + 
  geom_point(position=position_dodge(0.25))+
  geom_errorbar(aes(ymin=plb, ymax=pub), width=.2,
                position=position_dodge(0.25))
write_plot = p+labs(title="Write to Elected Official", x="Treatment Group", y = "Probability of Writing")+
  theme_classic() +
  scale_color_manual(values=c('black','gray')) +
  labs(col = "GPT3 vs Human")
write_plot

ggsave('write_plot.pdf', write_plot)
###################################################################